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Theoretical differences in the growth of structure offer the possibility that we might distinguish 
between modified gravity theories of dark energy and ACDM. A significant impediment to applying 
current and prospective large scale galaxy and weak lensing surveys to this problem is that, while the 
mildly nonlinear regime is important, there is a lack of numerical simulations of nonlinear growth in 
modified gravity theories. A major question exists as to whether existing analytical fits, created using 
simulations of standard gravity, can be confidently applied. In this paper we address this, presenting 
results of N-body simulations of a variety of models where gravity is altered including the Dvali, 
Gabadadze and Porrati model. We consider modifications that alter the Poisson equation and also 
consider the presence of anisotropic shear stress that alters how particles respond to the gravitational 
potential gradient. We establish how well analytical fits of the matter power spectrum by Peacock 
and Dodds and Smith et al. are able to predict the nonlinear growth found in the simulations from 
2 = 50 up to today, and also consider implications for the weak lensing convergence power spectrum. 
We find that the analytical fits provide good agreement with the simulations, being within la of 
the simulation results for cases with and without anisotropic stress and for scale-dependent and 
independent modifications of the Poisson equation. No strong preference for either analytical fit is 
found. 
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I. INTRODUCTION 

A diverse range of observations are showing consistent 
evidence for the acceleration of the universe's expansion, 
and the presence of dark energy, for example supernovae 
observations [U, 0, Q , cosmic microwave background tem- 
perature and polarization fluctuations 0, 0| , large 
scale structure surveys [1, Q and baryon acoustic oscilla- 
tions [loj . 

Interpretation of Einstein's cosmological constant as a 
vacuum energy requires the value to be fine-tuned to far 
smaller than any theoretical expectation, e.g-fllj). and 
has forced the exploration of alternative theoretical ex- 
planations. Since precision measurements of gravity only 
exist for scales < 10^'^ m (e.g.fl^l), there is freedom to 
posit modifications of gravity acting on larger scales, such 
as those in [l^. 

While cosmological observations of dark energy proper- 
ties have so far focused on measurements of the homoge- 
neous background density and astrophysical correlations 
in the linear regime, both theoretical and observational 
vistas are now opening up that require a good under- 
standing of the growth of structure in the mildly nonlin- 
ear regime. Theoretically, measuring the growth of struc- 
ture might enable modified gravity theories to be distin- 
guished from a standard cosmological scenario (ACDM) 
with a cosmological constant, A, and cold dark matter 
(CDM) [13, d, [ill . Observationally, the next genera- 
tion of precision experiments will include weak lensing 
surveys, with several proposed large scale weak lensing 
experiments being developed in the coming decade, e.g. 
DUNE [13, JDEM/SNAP [18]] and LSST [1^. 

Weak lensing is a potentially powerful probe of the late 
time evolution of the Universe, sensitive not only to the 
background expansion, but also able to give two-point 
and higher statistical correlations of the densityfield [13] , 
potentially in tomographic rcdshift slices 2^, 



Modified gravity models can introduce extrinsic 
anisotropic shear stresses (see e.g. that modify the 

relationship between the weak lensing potential and the 
matter over-density that might be detectable by contrast- 
ing weak lensing with other large scale structure obser- 
vations [H,!!!]. 

Many large scale structure statistics can be related to 
the underlying matter power spectrum, with nonlinear 
evolution at small scales. For standard general relativity, 
a typical approach is to use analytical fits based on N- 
body simulations of ACDM [H, ^ and CDM with dark 
energy with an equation of state, w,wCDM [23, [H, 
scenarios to apply the nonlinear correction to a linear 
power spectrum. Simulations of modified gravity models 
are for the most part lacking, however. With the ex- 
ception of [H, [m, analyses often proceed by applying 
the ACDM based analytical nonlinear fits to modified 
linear power spectra, e.g. [l^, [s^l- Recently an analyti- 
cal approach to estimating nonlinear growth in modified 
gravity, including those with anisotropic stress, [1^ was 
proposed and it was noted that there were currently no 
simulations against which to test the ansatz. 

In this work we directly address to what degree the 
nonlinear fits developed for standard gravity can be uti- 
lized for modified gravity theories, and establish whether 
the paucity of simulations in modified gravity theories 
can be excused. We first consider the applicability of 
standard gravity nonlinear fits to modified gravity theo- 
ries in which just the Poisson equation is modified, con- 
sidering the 5D gravit y fo rm considered by |34l | , comple- 
menting the work of [30l . [3l| who considered nonlinear 
growth when a Yukawa-like gravitational coupling is in- 
troduced [11] . We then address the impact of anisotropic 
stress on nonlinear growth to assess if existing analyti- 
cal nonlinear fits are adequate to model evolution in these 
scenarios. We consider the nonlinear growth in the Dvali, 
Gabadadze, and Porrati (DGP) 5D model [ll] and in toy 
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models that contrast the effects of anisotropic stress with 
those of a modified Poisson equation. 

We first estabhsh the framework for investigating mod- 
ified gravity theories in Sec. Ull and outhne the spe- 
cific models we consider with scale-independent and de- 
pendent modifications and the presence and absence of 
anisotropic shear. The details of our simulations and im- 
plementation, including the two standard analyitic fits 
are presented in Sec. IIIII The approach to weak lensing 
is discussed in Sec. IIVI The results showing dimension- 
less power spectra and the success of analytic fits are 
discussed in Sec. [V] An overview of the conclusions and 
implications of our results is then presented in Sec. IVII 



II. MODIFIED GRAVITY THEORIES 

We first outline the effect that the general modifica- 
tions to gravity wc study have on the perturbed Ein- 
stein's equations. Following the notation of [s^, in the 
conformal Newtonian (or longitudinal) gauge, the metric 
is written as 

ds^ = a(T)2 [-(1 + 2^jj)dT^ + (1 - 2(l3)dxjdx^ (1) 

where a is the expansion factor, r is the conformal time, 
X is the comoving coordinate (j=l,2.3 spatial directions) 
and (j) a-iid V' a-i'c the two gravitational metric perturba- 
tions. 

Einstein's equations relate the metric perturbations to 
fractional perturbations in density, Sg = 5ps/ps, peculiar 
velocity, W(s), and intrinsic shear as for a matter compo- 
nent "s", 



k^q^ + 3H{<}> + H^) = —QJ2^sSs, 



(2) 



k'i^ + Hi;) = ^^(l + z^,)r!,(zFi,(,),),(3) 



(t>-lp ^ — ^(l + Ws)l^sCr, +cro, (4) 

s 

where Ti, = a/ a, Q.{a) is the fractional energy density, and 
w{a) is the equation of state for the fiuid. We have in- 
troduced the function Q(k,a) as a modification in the re- 
lationship between the gravitational potentials and mat- 
ter density in the STq' equation. 121), and aQ{k,a) as an 
extrinsic anisotropic stress in addition to the intrinsic 
anisotropic stresses from the matter components (pre- 
dominantly radiation) in the equation for 5T^\ i ^ 3^ 
dH). For standard gravity Q = 1 and cto = 0. 
Equations ^ and ^ combine to give 



(6) 



where As is a gauge invariant density variable defined in 
the rest frame of the matter components f37| . 

Density and velocity perturbations evolve according to 
the perturbed fiuid equations which are unchanged by the 
gravitational modifications, 

S = -{l + w)iik^Vj -3^) ~micl~w)6, (7) 



iyij = — 



H{l-3w) 
k^a + k'^-ip, 



ik^Vi 



-k^S 



(8) 



where cj. is the sound speed for the fluid. 

We will consider a Universe dominated by pressureless 
matter, = = (Ts = 0, and scenarios in which ip ^ (f>, 
so that on subhorizon scales jfe^Vl ^ iSHt/il, and 

6 + n6 + k^ij « 0. (9) 

Using ([8]), we define the peculiar acceleration, g, 

1 d 



9j 



. dr 



(avj) = —ikjip. 



(10) 



Following the notation of [T^ , we relate the anisotropic 
stress to 4> through a function ij, 



(11) 



Q and r] here are equivalent to q and 7] in |38| . 

Making a subhorizon approximation, and Tiv/k <^ S, 
assuming v < 5, the modified Poisson equation and pe- 
culiar acceleration equations are 



37^2 

= ^Q^mS, 

9] = -ifcj(l + ?7)(^. 
while the matter perturbation equation is 
3>/2 

5 + m - —Q{1 + Ti)nm5 = 0. 



(12) 
(13) 



(14) 



Note that, we can describe the evolution of 6 in terms of 
the linear growth factor, D, with respect to some refer- 
ence scale, Oi, 5{k,a) = D{a)5(k,ai) where D is scale- 
independent for standard gravity, but could be scale- 
dependent if gravity is so modified. 

We can relate the Fourier space modification to a real 
space interaction in the form of a Green's function, 

0(r) = -Gpm{a)a^ J d\'Sir')f{r-r'), (15) 

g{r) = -V[(l + rKr))0(r)], (16) 

with /(r) = l/|r| recovering standard gravity. Using the 
convolution theorem wc find, 



Q{k, a) 



4n 



f{k,a). 



(17) 
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The effect of modified gravity in weak lensing statis- 
tics is described in [3^ where they show that the weak 
lensing distortion is dependent upon the sum of the two 
gravitational potentials, ^ = (j) + ip. As in [l^, we intro- 
duce the parameter £((5,77) to describe the deviation of 
the weak lensing potential from standard gravity 



E. ^1 
with S = 1 for standard gravity. 



(18) 
(19) 



A. 5D Gravity 

We consider a model, motivated by 5-dimcnsional grav- 
ity theories in which gravity is Newtonian on small scales 
but modified on scales larger than a characteristic scale 
rs 0, [13, Ell . This model is characterized by the form 



1 



(20) 



and 



Q(k,a) 



■ coa{krs) 




,(21) 



with r]{k, a) = 0. 

We are principally interested in the effect that mod- 
ifications to gravity could have on the transition from 
linear to nonlinear regime, typically occurring over co- 
moving scales 1 — 10 Mpc. For our analysis, therefore, 
we consider evolution for values of the parameter of 
20/i~i Mpc, lO/i^i Mpc, and bh'^ Mpc, which alters 
the behavior in the relevant scales. We do not consider 
here smaller values of the modification which would al- 
ter behavior in the wholly nonlinear regime. We leave 
it for future study to assess whether such changes are 
well modeled by analytical fits describing the properties 
of collapsed halos. 



B. DGP 

A physical model that serves as an excellent example of 
the effects of anisotropic shear is the Dvali, Gabadadze, 
and Porrati (DGP) model [3l that is based on 5D grav- 
ity, wherein at some large scale, Tc (comparable to the 
horizon scale), gravity is sensitive to the presence of an 
additional dimension. 

The extra dimension alters the 4D background evolu- 
tion to that described by the modified Friedmann equa- 
tion. 




, SttG , , , , 

TT-] +^Pm(a), (22) 



where H = W/a, with late time acceleration being trig- 
gered when the Universe's horizon ~ Tc- 

The modification also alters the growth of fluctuations 
in density and motion of matter. As well as a modifica- 
tion to the Poisson equation as discussed in Sec. Ill Al 
this model also results in an anisotropic shear such that 
the two potentials are given by [3, SI, S, IS] : 

ire 



2 



1 



1 + T-; 



1 

3/3 
1 

3"^ 



where 



/3 



(23) 
(24) 

(25) 



2H{a)rc - 1 

In contrast to Sec. Ill Al this gives a scale-independent 
modification to the Poisson equation. 



Qia) ^ 1-^, 

and nonnegligiblc anisotropic stress, 

2 



via) 



3/3-1 



(26) 



(27) 



and E = 1. 

For our analysis, with a background cosmology with 
Hubble constant, i?o,= 70 km s^^ Mpc^^, fractional 
matter density, flrn,= 0.3, and consistent with the ob- 
servational constraints found in [3| (EH), fc = 6.1 Gpc. 



C. Twin toy models 

Finally, we consider a set of twin models that provide 
a simple way to further explore the effects of anisotropic 
stress on nonlinear growth. We consider two different 
modifications that both yield the same form for the weak 
lensing potential, with E = l-|-Eoa, such that they reduce 
to standard gravity at early times and become modified 
at late times. This form of E is equivalent to model 
GDEl of [H. The twin models ("TM") we study have 
two contrasting, simple forms in terms of Q and rj : 



and 



TM 1 : Q = 1, r; = 2Soa, 



TM 2 : Q = 1 + Ena, 7? = 0. 



(28) 



(29) 



In TMl, the Poisson equation is the same as for stan- 
dard gravity; however, the peculiar acceleration of the 
matter particles responding to the gradient of the po- 
tential is affected by the anisotropic stress. In TM2, in 
contrast, the peculiar acceleration is the same as for stan- 
dard gravity but the gravitational potential at late times 
has a different relation to matter over/under densities. 
We consider values of Eq = ±0.008, ±0.016 consistent 
with 1 and 2a Fishcr-matrix constraints for a prospec- 
tive DUNE- like weak lensing survey p^ . 
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III. N-BODY SIMULATIONS 

To obtain fully nonlinear results in each of the models, 
we obtain N-body simulations via a particle mesh(PM) 
code, taking as an initial form the code of [iH . For scale- 
independent modifications we make simple modifications 
to the code, described in Sec. IIII Al For scale-dependent 
modifications we have to alter the potential and motion 
calculations as described in Sec. IIIIBl 



^1 Cell(Ij- l,J+l) 



4^< — > 



Cell (I+1,J) 




A. Standard Gravity and Scale-independent 
Modifications 



The PM code is reviewed in detail in [45 but we pro- 
vide some highlights in order to set the framework for 
discussing the modifications we make to the code. 

PM codes operate by defining a simulation area as a 
box of size L on a side, assuming it is closed so that we 
have periodic boundary conditions, subdividing it into a 
mesh or grid of cells (of size L/N on a side), and defin- 
ing all quantities on that mesh. The simulation is then 
initialized at some early redshift (z^) and Np particles 
are placed according to model-dependent power spectra 
fits provided with the code (based on the cosmological pa- 
rameters: the scalar spectral index n^; the amplitude of 
fluctuations in 8h~^ Mpc, (Tg; the fractional density from 
curvature, fix, baryons, fib and cold dark matter, ^cdm', 
and the Hubble constant Hq = lOOh km s~^ Mpc~^). 
The evolution is then carried out by advancing time in 
equal steps of the expansion factor, Ustep- At each step 
in expansion factor the code determines a density in each 
cell, uses that density to compute the potential </) in each 
cell, and finally moves particles according to the gradient 
of the potential. 



1. Defining the Density 

Defining the density can be done in a variety of ways; 
the code uses the cloud-in-cell scheme depicted in Fig. [T] 
wherein a particle is taken to be a cube with dimensions 
equal to that of the cells and with a corner positioned 
at the location of the particle. The particle contributes 
to each cell it extends into a mass equal to the particle's 
total mass weighted by the fraction of the particle's vol- 
ume in the cell under consideration. Once the mass in 
each cell is determined it is effectively smeared over the 
entire cell. 



2. Obtaining the Potential 

For standard gravity, the code us es (fT2l) with Q = 1, 
with the dimcnsionlcss variables of [45| . f = r/xo and 



FIG. 1: A two dimensional description of cloud in cell density 
assignment, (a) The definition of the variables in relation 
to the particle's actual position. The particle is the black 
dot, but it is extended to be a square particle denoted by the 
dotted lines, thus it lies in four cells. The sides of the cells and 
the size of the particle square are L = Dl + Tl = D2 + T2. 
(b) The resultant mass distribution in each cell. Note that 
the mass is not retained in the original particle's area, but 
rather smeared over the cell it occupies. 



(/) = (/)/{xoHo)'^ and writing S = p{x,a)/p{a) — 1, 



3 ^rn.Q ^ 

2 a ' 



(30) 



To evaluate ([50]) , we use the discretized Poisson equation 
over cells, n = 0, iV — 1. In one dimension, the discrete 
Laplacian is given by 



Defining the discrete Fourier transform, 

<Pk — ^„=o Vne ' , 

the discretized Poisson equation is 



(31) 
(32) 

(33) 



Generalizing to three dimensions one obtains the '7-point 
crest template', 

V^0ij-,fc ~ + 4>i-lJ,k + 4>i,j + l,k + 4>i,j-l,k 

+ 4>i,jM+i + 4>i,j.k-i - &4>i,j,k, (34) 







- 1 


X 2 


cos 




. \N ) 





with 

where Gk is given by 

' \ ( 2'nk. 



(35) 



Gk 



N 



N 



2Trk^ 



-3 
(36) 

Combining ((35|) and ((30)) the Poisson equation used in 
the code is, 



7 3 flm,0 r 



(37) 
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The code calculates ^(r), Fourier transforms to (5(k). di- 
vides by Gfe and then transforms back to real space to 
obtain 4>{i,j, k). 

In the case of scale-indepcdent modifications ([ST]) is 
purely modified by 

GkMk,ci) = ^. (38) 



3. Advancing the Particles 

Once we have the potential (j) we advance the particles 
according to (|16[) . In standard gravity, component wise 
on the grid we have only to compute 

9x = -{4>i+i.].k - 0i-ij,fc)/2 

9y = -{(f'ij+lM ~ 0ij-l,fe)/2 

In the presence of anisotropic stress modifications, 

gj,ait = [l + '?(a)]3j- (39) 

B. Scale-dependent Modifications 



The final 3-dimensional radius, f, is computed trivially 
as 

r'=j2rl (41) 

i=l 

There remains one final subtlety in computing the ra- 
dius. Since f[l,l,l] = 0, division by r requires us to 
make a change to avoid infinities. To avoid these sin- 
gularities we take the standard approach of 'softening' 
f (e. g. p6|). that is adding a small non-zero term to 
all the values of r used in operations that would give 
a singularity. For instance if we consider g{f) = 1/f we 
instead compute g{f) = (f^ +e^)~^/^. Note that, for con- 
sistency, all values of r in the division are softened, not 
only the actual one that gives a singularity (f[l,l,l]). 
Further, note in the case of well defined modifications, 
e.g. g{r) = e~^/r the exponent need not be softened, so 
that we compute e~''(f^ -|- e^)~^/^. 

C. Obtaining Analytic Spectra 

We compare the nonlinear spectra from simulations 
to predicted spectra from analytical mappings of linear 
power spectra using the Peacock and Dodds (PD) fit [1^ 
and the Smith et al. fit (SP) [H]. 



In order to incorporate the scale-dependent modifica- 
tions to gravity we follow the convolution approach in 
([T5|l . To do this we multiply by /(k, a) at each step in a 
rather than 1/Gk- 



1. Defining the Radius r for g{r) 

In scale-dependent theories, by definition, we now con- 
volve with functions involving the actual scale r, and we 
must therefore define explicitly a radius on the grid. The 
mass is smeared over the entire cell it lies in, so that the 
distances simply become those between cells. The Fourier 
transforms involve periodic boundary conditions, so we 
define the radius for one origin at (0, 0, 0), and wrap the 
radius around the grid. Since the code uses the dimen- 
sionless radii to compute the function we have called /(f), 
the cell indices can be used to construct the radius and 
we define fi to be the index of the relevant cell in the i*'' 
direction (i=l,3). 

The periodic boundary conditions require a change to 
the basic prescription presented above, namely to include 
the periodic boundary condition we must set up a ID 
radius of the form [0, 1, 2, N/2 - 1, N/2, N/2 - 
1, 2, 1] where N is the number of cells making up the 
grid in a dimension. Thus, the radii in the i-th dimension 
can be defined as 

~ / \ f n< N,/2 

'^'(")-|iV.-r. n>N:/2 ^^0) 



1. Analytical Linear Spectra 

We evolve a linear ACDM power spectrum obtained 
with CAMB [131 (that includes effects from baryon pho- 
ton coupling at early times) forward in time using the 
modified equation for the growth of the over-density (|14p . 
We start at an epoch, we choose Zi = 50, at which the 
modification scale is large compared to the physical hori- 
zon, so that standard gravity is effectively recovered on 
the relevant scales, and evolve the density perturbations 
through the modified gravity era to today. In Fig. [2] we 
show the linear power spectra for the models discussed 
in the paper. 

2. Analytical Non-Linear Fits 

We briefly review the physical ingredients of the PD 
[25t and SP [1^ analytical fits against which we compare 
the simulations. 

The PD fit is based on the assumption of stable cluster- 
ing [i^, the hypothesis that the correlation function on 
scales smaller than those of virialized structures decouple 
from the expansion. The fit utilizes a linear to nonlinear 
mapping proposed by Hamilton et al. (HKML) [4§| 

kL = [l + A%L{kNL)r'^^kNL. (42) 

derived from the spherical collapse model. Peacock and 
Dodds generalized the HKML method to estimate the 
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TM1 0.016 TM2 0.016 

TM1 -0.016 TM2 -0.016 

5-D Gravity Model r^=20 h ' Mpc 
5-D Gravity IVIodel r_,=10 h"' IVIpc 
5-D Gravity IVIodel r^=5 h ' Mpc 
DGP r =6.1 Gpc 



k(Mpc"') 



FIG. 2: The ratio of the hnear power spectrum in the modified 
theories to that for standard gravity for the models discussed 
in Sec. [Ill the 5-D gravity model of Uzan and Bernadeau 
(solid line), TMl (dotted line), TM2(dashed line) and DGP 
(dotted-dashed line). 



resulting nonlinear power spectrum through a universal 
scaling relation, /jvl, 



fNL{x) = X 



1 + Bl3x + [A: 



,10/3 



l + {[Ax]"g^a)/[Vx^/^])P 



(43) 
(44) 



where g = D{a)/a. The fitting function /at^ tends to 
fNL{x) = a; in the linear limit, x ^ 1, and fNL{x) = 
V g~^{Vlm)x^/'^ in the small scale, stable clustering limit, 
X ^ 1. There are five free parameters fit from N-body 
simulations in standard gravity as functions of the linear 
spectral index rie// = dlnP^ / d\Tik{k = A;i/2): A and a 
parameterize the power law in the quasilinear, large scale 
regime, V parameterizes the amplitude of the fNL{x) in 
the stable clustering limit, B describes the second order 
deviation from linear growth and P softens the transition 
between the linear and fully virialized regimes. The cos- 
mological model only enters into the fit through g, con- 
sistent with the Zel'dovich approximation in which the 
final positions of particles are obtained by extrapolating 
their initial comoving displacements, q, using the linear 
growth factor, x{a,t) = a{t) [q + D{a)Vijj{q)]. 

The quality of the PD fit is founded on the broad appli- 
cability of the Zel'dovich approximation. However, with 
a scale-dependent modification of gravity, or the intro- 
duction of a difference between and ip it is not clear a 
priori how well the Zel'dovich approximation will apply, 
and if applicable, whether the numerical values of the 
coefficients will remain the same as those for standard 
gravity. That is, with scale-dependent modifications the 
possibility for shell crossings arises which causes a break 



down of the Zel'dovich approximation. 

Looking at the functional form of the fit, in particu- 
lar, we can consider three regimes to make predictions, 
namely, the large and small scale limits and a transition 
regime. Large scales which remain linear or quasilinear 
should be well described by the existing fit. On these 
scales the Zel'dovich approximation should hold and us- 
ing a linear growth factor for g is acceptable. Similarly a 
and f3 might be expected to adapt to the changed input 
power via their spectral index dependence, since in linear 
scales essentially all the information is contained in the 
amplitude and spectral index of the power spectrum. 

The mildly nonlinear or transition regime, where we 
directly compare results, is particularly of interest in ap- 
plying the fits. Scale-dependent modifications introduce 
an extra degree of freedom to growth in the model, a scale 
dependency that could also affect the shape and scale of 
the smoothing function interpolating between the linear 
and nonlinear asymptotic behaviors, essentially requiring 
corrections to /?. For example, a scale- or time-dependent 
modification to Poisson's equation could alter the critical 
over-density required for nonlinear collapse, thus alter- 
ing the details of the transition from linear to nonlinear 
regimes. 

Small scale modifications to gravity, which we do not 
consider here, could well lead to alterations in the corre- 
lation function of the collapsed structures, in particular, 
changes to the value of the virialized normailzation V. 
One might expect the application of the linear growth 
factor in the fit to be less effective even if including the 
linear scale dependency g{a) g{k,a). Relevant to our 
analysis is the fact that the stable clustering approxima- 
tion does not account for merging and accretion of halos 
and hence does not address how modifications to gravity 
may alter these physical processes. We discuss this in the 
context of the SP fit below. 

The SP fit arises from a different approach based on the 
"halo model" [10, HH in which the continuous accretion 
of matter and merging of halos is accounted for, deviating 
away from the stable clustering approximation. In this 
scenario, the power spectrum of matter A^^ = Ag -I- 
is described on large scales by the correlations between 
different halos represented by a quasilinear term, AQ(fc), 
and on small scales by a halo term, A|^(fc), that accounts 
for power from the self-correlation of halos. In the fit, 
the two terms arc phenomcnologically selected functions 
oi y = k/k^, where the scale k^{a) becomes nonlinear at 
scale factor a(t). 



A|(fc) = Ai(fc) 
Alik) = 



(1 + Ai(fc)) 



/3„ 



any 



3/i(n„ 



l + a„A2(fc) 



cxp 



y , y 



(45) 



l + 5„2;/2(0„) + (c„/3(a„)y)3-^" 

1 

1 + fj-ny'^ + Vny'"^ ' 
k„ is determined by the standard error of the linear den 



(46) 
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sity field, 

<T{k-\a) 
<t{R, a) 



1, 



J Al{k,a)cxp{~k^R^)d\nk. 



(47) 
(48) 



The eight coefficients {a„, /3„, 7„, i^„, a„, c„}, fit 
with spectral index-dependent functions, and three $7^- 
dependent functions, /i, /2, and /a, are empirically 
matched to standard gravity simulations. 

In the large scale limit, the quasilinear term dominates 
and the use of the spectral index dependent an and /3„ 
suggest the fit will adapt well to a modification on linear 
scales, such as those considered here. 

On small scales, just as in the PD case, there are is- 
sues with the numerical fitting functions in the halo self- 
correlation term; the correlation coefficient in a virialized 
halo could be modified for the various modified gravity 
scenarios. The functions fi, /2, and /s, which in standard 
gravity are purely functions of 17™, would be expected to 
alter to account for the modification; this in turn could 
well be expected to change a„, 6„, c„. 

Also as in the PD case, the interpolation from linear 
to nonlinear regimes, from large to small scales, could be 
altered as the modifications could alter the critical over- 
density required for nonlinear collapse. In particular the 
form of a„ and to some extent /3„ may be expected to 
require changes as these serve to determine the relative 
importance of the halo- halo and self-correlation terms. 

It is in light of these considerations that we study 
whether these analytic fits can readily describe modi- 
fied gravity scenarios, with scale-dependent or indepen- 
dent modifications to the Poisson equation, and/or scale- 
independent anisotropic shear. 



IV. OBTAINING WEAK LENSING SPECTRA 

Modified gravity theories can impact weak lensing con- 
vergence power spectrum in addition to the matter power 
spectrum, thus we study the impact of our models on 
both. In standard gravity, the power spectrum of the 
convergence is given by 



40"* 4 Jo a^X 



where Ps is the matter power spectrum and g{x) is a- 
weighting function that can be related to the comoving 
distance x a-nd the distribution of background or source 
galaxies, Ws(x) 



X 



X 



(50) 



We assume a a simple delta function distribution of 
sources at = 1, so. 



9{x) = 2x 



Xs-X 



(51) 



The convergence power spectrum is then 



4c3 a'^H{a)x{af 



Xia) 



with 



W{x,Xs) = X 



Xs-X 

Xs 



da (52) 



(53) 



Where we have used the fact that the comoving distance, 
X, is equal to the (comoving) angular diameter distance 
for a fiat Universe so that 



xia) 



cAa' 



(54) 



Gravitational modifications Q ^ 1 and/or rj ^ 0, will act 
to modify Ps(k, a). In addition, Q ^ 1 and/or ?/ ^ will 
modify how the convergence spectrum is related to Pg 
(see for example (isj). resulting in 



4c3 



' wHx,XsKi + v/2 ) 

a^H{a)x{aY 



2n2 



-Ps 



xia) 



da, (55) 



where the evolution of Ps is also affected by Q and 77. 

Ps can be obtained from either the PM simulations or 
the analytic fits described in Sec. IIII CI To actually eval- 
uate the integral, we discretize it, binning by expansion 
factor. 

When considering the N-body code derived Ps , wc have 
to account for the fact that the simulation only probes 
a range of fc, yet for any given I, k = I /xia) can can lie 
outside this range at some redshift, Zg > z > 0. For the 
I range we consider, the range of k needed is virtually all 
given by the N-body simulation. Outside this range, on 
large scales the power spectrum is well approximated by 
the linear spectrum; at smaller scales we find that the 
analytic SP and PD predictions for the modified gravity 
spectra are within the Icr errors at the edges of the range 
of k provided by the simulations, so we pad the simu- 
lated spectra with the nonlinear analytic fits to modified 
gravity linear power spectrum. 



V. RESULTS 

A. Parameters for Simulations 

For the PM code parameters wc take N = 256, Np ~ 
128, L = lOO/i-i Mpc, e = 0.1, Zi = 50, astep = 0.002 
and for our cosmological model we take = 1, erg = 
1, Qk = 0, rirn = 0.3, nb = 0.026, ricdm = 0.274, 
and h = 0.7. The resulting simulations measure scales 
0-1 ^ ^ ^ 1 Mpc^^. The specific choice of initial redshift 
is not important other than to ensure that it is early 
enough that nonlinear corrections are negligible. 
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The box size and number of cells play into spatial res- 
olution of the simulation, and are chosen to allow us to 
effectively probe the decade of k in which the mildly 
nonlinear effects manifest themselves and from which 
we can extract a reasonable weak lensing spectrum for 
I ^ 200 — 1000, a range relevant to upcoming experi- 
ments. The number of particles are chosen to ensure a 
sufficient particle resolution for the box size and number 
of cells used. 

The initial positions of the particles at Zi , are assigned 
by means of a random number generator consistent with 
the initial power spectrum. Depending on the seed used 
to initialize the random number generator, the resultant 
spectra may agree well with standard ACDM analytic 
fits with the same parameters or might over- or under- 
produce power, even in the original unaltered code of [i^ . 
We therefore run the simulations with 24 random seeds 
to get a good sample size and a more robust average. In 
order to weight the behavior of each simulation equally, 
we consider the modifications in terms of the "average of 
the ratios" of the modified power spectrum to the stan- 
dard gravity spectrum for the same seed, rather than the 
"ratio of the averages" that would preferentially weight 
those simulations that over-produce power. 

For scale-dependent modified gravity, we find the Nu- 
merical Recipes routine (52j for the Fourier transform, 
though slightly more time consuming, is more stable than 
the one provided in the original code. In the case of scale- 
independent modifications and standard gravity, both al- 
gorithms produce identical results. The softening param- 
eter value used for the scale-dependent modification is 
much smaller than the smallest separation in the code 
and provides agreement with standard gravity from an- 
alytic predictions and standard gravity simulations with 
the code at least at the level or better than the unmodi- 
fied Klypin code. 

In Fig. [3] we show the results of the 24 simulations of 
standard gravity against the SP and PD fits, in order to 
demonstrate the fiducial model against which the modi- 
fied gravity simulations are compared. The simulations 
are consistent with the analytical fits in the range 0.1 
Mpc~^ ^ k < 1 Mpc~^. A conservative estimate for the 
largest k at which we can believe the simulation results 
are reasonable is k Nyquist [lOl, which for our simula- 
tions is 1.4 Mpc~^. We consider the simulations to be 
valid only up to kj^yquist/'^, rather than up to k^yquist as 
this more conservative limit represents a regime in which 
standard gravity simulations and fits agree to within 1.5 
times the standard error in the simulation, in comparison 
to 10 (for the PD fit) and 13 (for the SP fit) times the 
standard error at kjsiyquist- 

For the model parameterizations we consider, we find 
that the linear scales used to generate the nonlinear k in 
the range 0.1 — 1.4 Mpc~^ lie in the range k ~ 0.07 — 0.5 
Mpc-^ 
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FIG. 3: Dimensionless matter power spectrum, A^(fc) = 
fc'^Pi-(fc)/27r^, for standard gravity. The full line and errors 
bars show the average power spectrum and standard devia- 
tion for 24 simulations. The vertical dotted line represents 
kNyquiat/2, which is a conservative estimate for the largest k 
at which we can believe the simulation results as in [s^l ■ The 
PD (dot-dashed line) and SP (dashed line) analytical fits are 
also shown. 

B. Simulation and Analytic Fit Results 

1. 5D Gravity Model 

The ratio of the dimensionless power spectrum today 
for the 5D gravity model discussed in Sec. Ill Al to stan- 
dard gravity, is shown in Fig. 31 We find that the simu- 
lations are consistent with the PD [2^ fit at the la level. 
This is consistent with the results of [l^l for a Yukawa 
type modification (that, like the modification we consider 
here, is a scale-dependent modification). The SP fits are 
slightly less consistent with the numerical predictions, 
however, still lie within 1 cr of the simulation mean. We, 
therefore, find no statistical basis for preferring PD over 
the SP [2| fit ofPsiz = 0). 

To consider the suitability of the analytic fitting func- 
tions when applied to weak lensing, it is insufficient to 
purely consider their agreement with predictions today; 
the entire evolution must be tracked between the redshift 
of the lensed source and today, as weak lensing integrates 
P{k, a) over the expansion factor a, c.f. ([55)1 . We, there- 
fore, track the redshift history of the nonlinear evolution, 
and the comparison with the analytical fits, as shown in 
Fig. O Wc find both fits lie within Icr though after 
a ~ 0.75 the SP results are just encompassed by the la 
errors. 

The ratios of the modified gravity weak lensing spec- 
tra to those of standard gravity are well recovered by 
the PD and SP fits, as shown in Fig. [S] The ratios 
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FIG. 4; Ratios of the z = dimensionless matter power 
spectrum in the modified gravity model to that for standard 
gravity, for the 5D gravity model described in Sec. Ill Al for 
rs = 20h~^ Mpc (top, blue), lOh^^ Mpc (middle, green) and 
5h~^ Mpc (bottom, red). The full line and errors bars show 
the average of the ratios and standard deviation for 24 simu- 
lations. The vertical dotted line represents kNyquist /2, which 
is a conservative estimate for the largest k at which we can 
believe the simulation results as in [30]. The PD (dot-dashed 
line) and SP (dashed line) analytical fits agree with simula- 
tions to within la for each r^, in the region of interest, k = 0.1 
to 1 Mpc"\ 



of the weak lensing convergence spectra are slightly less 
sensitive to the exact form of the modification than the 
matter power spectra, for two reasons. First, the inte- 
gral in (j55p is mostly weighted towards integrand values 
at early times when the analytical fits are in very strong 
agreement with the simulations. Thus, for instance, the 
late-time transition of the SP fit to the outer regions of 
the 1-sigma level is not so significant to the convergence 
power as it is to the final matter power spectrum. Sec- 
ondly, we "pad" the spectrum at k values outside the 
simulated range with the analytical fit, in order to eval- 
uate / = kx{a) in (|55p . This is mitigated (as discussed in 
Sec. IIV[) by the similarity of the fits and the code spectra 
at the edges of our range of k and the fact that the con- 
tribution from padded k values is small in comparison 
to those drawn from the simulated range: for I — 200, 
Pk is padded with the nonlinear analytical spectrum at 
a > 0.955, which corresponds to 1.3% of for stan- 
dard gravity; for I ~ 1000 the padding is required for 
0.8 < a < 1.0 which contributes to 14% of the value of 

Pn- 



FIG. 5: The ratios of the dimensionless matter power spec- 
trum in modified to standard gravity, A'^n{k) / A'l-i.j^{k) as a 
function of redshift 50 < z < for A; = 0.53 Mpc"^ The 
color coding and lines styles are as in Fig. |4l The dotted 
lines show the ratios of the associated linear spectra. Note 
that the evolution is well tracked by the analytical fits, with 
both lying within la for the simulations. At late times the 
SP fit drifts to around, or just over, the la error. 



2. DGP 

The effects of nonlinear growth in DGP models are of 
great interest in establishing observational distinctions 
between this model and standard ACDM at cosmologi- 
cal scales, in jisj the nonlinear power spectrum was esti- 
mated using the Smith et al. analytical fit, while in (ssj 
an analytical ansatz is applied. Both the DGP model and 
the model in Sec. Ill Al are motivated by 5D modifications 
to gravity. The difference between DGP and that model 
is that DGP not only modifies the Poisson equation but 
also the peculiar acceleration through the presence of an 
anisotropic stress. 

For the arguably more complex DGP model, the SP 
and PD fits are both still in good agreement with the 
N-body simulations at a = 1, at the la level over the 
simulated scales, as shown in Fig. [71 This is also true 
over the course of the evolution as the modification from 
ACDM switches on, as shown in Fig. [5| 

Note that we do not provide a weak lensing analysis 
in this model; as due to the change in H [and hence in 
x(a)] evaluating k — l/xio) consistently results in a need 
for much smaller scales, i.e. k > 4.6 Mpc~^ for the range 
of Vs we have considered. We thus restrict our discussion 
of DGP to matter power spectra and their evolution. 

Even though rc is chosen to be in close agreement 
with the background evolution of our fiducial cosmo- 
logical model, and has essentially degenerate evolution 
at early times, the DGP model shows marked deviation 
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FIG. 6: The ratio of the weak lensing dimensionless convergence power spectrum, = l'^ Pk{1)/'2'jt, for a S function lensing 

source at Zs = 1, as a function of muhipole, I, for the 5D gravity model in Sec. Ill Al to that in standard gravity in comparison to 
the SP fit (left-hand panel) and PD fit (right-hand panel). The points and errors are the average and standard deviation of the 
ratios the 24 simulations. The predicted spectra from the analytical fits (full lines) are wholly consistent with the simulations 
for all 3 modified gravity models with = 20h^^ Mpc (top, blue), lOh^^ Mpc (middle, green) and 5h^^ Mpc (bottom, red). 



from standard gravity at late times. Wc note that the 
suppression of the nonlinear power spectrum shown with 
respect to standard gravity for the PD and SP fits and 
N-body simulations is qualitatively similar to that shown 
with the ansatz of [33j , although we leave a quantitative 
assessment of the ansatz to future work. 



3. Twin Toy Models 

In order to investigate the abilities of the two analyt- 
ical fits to predict nonlinear behavior in the two types 
of modifications, we consider a set of twin toy mod- 
els, described in Sec. Ill CI TMl has a modified Pois- 
son equation {Q = 1 + Soa, '7 = 1} while TM2 has 
anisotropic stress {Q = 1,?/ = 2Y,Qa}. Both models 
have the same form of relationship of the weak lensing 
potential to the over-density, characterized by the func- 
tion E(fc,a) = Q{1 + ry/2). As shown in Fig. [TTl de- 
spite the degenerate background evolutions, the different 
modifications in each model lead to different linear scale- 
independent growth factors. For both models, the SP and 
PD analytical fits track both the scale-dependent behav- 
ior and time evolution of nonlinearities in both types of 
scenario, as shown in Figs. P and [TUl The weak lensing 
correlations for PD and SP fits are virtually identical for 
each model so we only show the results for SP fits in Fig. 
nil the difference between the simulations and analytical 
fits is negligible for both models. 



C. Discussion 

The nonlinear fits of Peacock and Dodds and Smith et 
al. have been shown to work across broad cosmological 
models with standard gravity, with different fractional 
mass densities, curvature, and initial power spectrum 
spectral indices. The utility of these fits derives from 
the wide applicability of the Zel'dovich approximation. 
In both fits there is the conjecture that the statistics of 
the gravitational clustering obey a similarity transform 
PNL(k/a) = P{k/kpfL) for which no proof is given, but 
instead is experimentally shown to be robust for a variety 
of cosmological models by simulation. In this paper we 
assess whether such a similarity transform similarly ex- 
ists in modified gravity theories, and moreover that the 
existing quantitative values for the fit coefhcients can be 
used. This is not necessarily the case a priori. 

To test the fits we have performed nonlinear simula- 
tions of models in which modifications to Poisson and 
the peculiar acceleration occur exactly in this mildly non- 
linear, transition regime. We have found that both the 
SP and PD analytical fits give reasonably good agree- 
ment with the simulations, in spite of the scale- and 
time-dependent modifications. This implies that appli- 
cability of the Zel'dovich approximation extends to sce- 
narios in which anisotropic stress and even those with 
scale-dependent modifications to gravity are present in 
the mildly nonlinear regime. The modifications, there- 
fore, are well described by the fits through their impact 
on the linear growth factor, (7, and the spectral index de- 
pendency of the fitting functions. It appears that scale- 
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FIG. 7: Ratios of the matter power spectrum in the DGP 
model with —6.1 Gpc~^ to that in standard gravity; botli 
models have Hq = 70 km s"^ Mpc"^ and Q.^ = 0.3. The fuU 
line is the average of the 24 realizations and errors represent 
the standard deviation of the simulations. The SP (dashed 
line) and PD (dotted-dashed line) analytic fits are in good 
agreement over the scales measured by the simulation, = 0.1 
to 1 Mpc~^. The linear power spectrum ratio is shown by the 
dotted line. 

dependent modifications in the mildly nonlinear regime 
do not require significant modification of the numerical 
coefficients in the fitting functions. Since our simulations 
focus on the ability of the fits to accurately match the 
transition from linear to nonlinear regimes, they do not 
investigate if modifications on small scales, in which the 
subhalo correlations are key, are well described by the 
fits, for example, if in (j20p were significantly smaller, 
e.g. less than 1 Mpc. This is an area of interest for fur- 
ther analysis, especially in recently discussed theories in 
which galactic scale modifications could be present (e.g. 



VI. CONCLUSIONS 

The use of complementary cosmological observations 
to probe the properties of dark energy has proved ex- 
tremely powerful. Observations sensitive to the back- 
ground evolution, e.g. supernovac, or the wholly lin- 
ear regime, e.g. the cosmic microwave background, have 
been the major observational tools to constrain dark en- 
ergy to date. There is now significant interest, however, 
in applying a broader range of observations including 
those sensitive to large scale structure including large 
scale galaxy surveys, such as the Sloan Digital Sky Sur- 
vey, and current and prospective weak lensing surveys. 
For each of these, in order to make precise inferences 



FIG. 8: The evolution of the ratio of the DGP matter power 
spectrum to standard gravity for k — 0.53 Mpc~^ as a func- 
tion of scale factor, a. The full line is the average of the 24 
realizations and errors represent the standard deviation of the 
simulations. The SP (dotted line) and PD (dotted-dashed 
line) analytic fits are good at predicting the transition and 
development of nonlinear growth at all epochs. 



about dark energy, theoretical systematic errors about 
the modeling of nonlinear corrections must be addressed. 

This work considers the ability of the commonly used 
nonlinear analytical fits of Peacock and Dodds [2^ and 
Smith et al. |26j to predict nonlinear growth in a va- 
riety of theories beyond standard gravity. We consider 
models in which the Poisson equation is modified, based 
on 5D gravity, [13, SO, EH and also those in which pecu- 
liar acceleration response to the gravitational potential 
is altered, including the DGP model [Tsj . 

We find that the two fitting functions provide robust 
predictions for theories with both types of modification, 
in terms of accurately predicting the matter power spec- 
trum today, and also, vitally for calculating the weak 
lensing convergence spectrum, they predict the develop- 
ment of nonlinearities over time. Both consistently give 
predictions within Icr of 24 simulated N-body realiza- 
tions of the theory. Our results imply that the similarity 
conjecture for mapping linear to nonlinear power empiri- 
cally found to be satisfied in standard gravity simulations 
is also applicable to scenarios in which gravity has scale- 
and time-dependent modifications. This suggests that 
the spectral index dependence of the fitting function and 
the linear growth factor effectively describe alterations 
in the nonlinear collapse due to scale-dependent modi- 
fications to gravity and anisotropic stress at the scales 
studied in the models here. 

We conclude that current analytic fits using the linear 
power spectrum in modified gravity theories can be used 
to accurately predict the nonlinear growth in theories 
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FIG. 9: The ratios of matter power spectra at a = 1 for modified gravity to standard gravity in the TMl(left panel) and TM2 
(right panel) models for Eo — —0.016 (dark blue, bottom), -0.008 (red) ,0.008 (green) and 0.016 (light blue,top) as a function 
of scale, k. As in earlier figures, the full line represents the average of the 24 simulations, error bars represent one standard 
deviation, and fcjvi/<jiiist/2 is indicated by the vertical dotted line. The predictions of the SP (dashed line) and PD (dotted- 
dashed line) fits are nearly identical, and are in excellent agreement with the simulations for both the weaker modifications 
with Eo = ±0.008 and the strong ones with Eo = ±0.016. The linear power spectra, showing the differences in linear growth 
factor arising from the modifications are shown by the dotted lines. 
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FIG. 10: The evolution of the power spectrum over time for TMl (left panel), and TM2 (right panel). Throughout the entire 
simulation the fits track the simulation results results extremely well. The color coding and line styles are the same as in Fig. 
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discussions in the course of this work. The work of IL 
and RB is supported by the National Science Foundation 
under Grants No. AST-0607018 and No. PHY-0555216. 



FIG. 11: The ratios of modified convergence power to stan- 
dard convergence power in the twin models TMl (full trian- 
gles) and TM2 (empty triangles) for Eq = 0.016 (<, blue), 
-0.008 (A,red) ,0.008 (V, green), and 0.016 (>, light blue) 
shown against the predicted spectrum using the SP fit (full 
line) , as the predictions of SP and PD are virtually identical. 
As is to be expected, given the strong agreement between the 
fits and simulations of the matter power spectrum, the weak 
lensing spectra from the simulations are predicted well by the 
analytical fits. 



with scale-independent or dependent modifications, and 
in those with or without anisotropic stress in the mildly 
nonlinear regime. We find no statistical evidence for a 
preference, on the basis of overall performance, for one 
analytical fit over the other. 

Many modified gravity models, for example, DGP and 
/(i?) models, exhibit gravitational modifications on sub- 
halo scales. Whether such modifications are well de- 
scribed by the halo term in the SP fit or the stable clus- 
tering approximation in the PD fit necessitates smaller 
scale simulations in the substantially nonlinear regime, 
which lies outside the scope of this paper. 

We have limited our investigation of anisotropic stress 
to scale-independent modifications, and indeed further 
work is warranted in investigating whether the conclu- 
sions found for those are applicable to scale-dependent 
anisotropic stress, as found in f{R) theories. It will also 
be interesting to investigate the agreement between sim- 
ulations and the recently proposed nonlinear ansatz [s^ 
for modified gravity models. 
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